Recent advances of 3D seismic survey led to the development of the data processing technique for 3D seismic data. In particular, the use of 3D reverse-time migration to study subsurface structures has attracted much attention (Abdelkhalek et al., 2009; Araya-Polo et al., 2009; Kim et al., 2011; Yoon et al., 2003). As the reverse-time migration method became more popular, a more accurate velocity model was needed. As a result, recent research in data processing areas has focused on accurate velocity model building, and full waveform inversion is one of the available techniques. Recently, research on 3D full waveform inversion is carried out by many geophysicists (Ben-Hadj-Ali et al., 2009; Plessix, 2009; Pyun et al., 2011b; Son et al., 2014). Under these circumstances, a practical problem has arisen in that the reverse-time migration requires enormous computational costs to verify the usefulness of inversion results. Therefore, a cost-effective migration technique is needed to verify the waveform inversion results. Although Kirchhoff migration is not as accurate as reverse-time migration, it is efficient enough to verify 3D inversion results.
To carry out Kirchhoff migration, underground traveltime information is necessary, wherein the traveltime can be calculated by many algorithms such as various ray tracing methods (Coultrip, 1993), Eikonal solvers (Vidale, 1988; Vidale, 1990) and wave-equation-based algorithms (Shin et al., 2002; Shin et al., 2003; Qin et al., 2005), etc. Although ray tracing methods or Eikonal solvers are more efficient than wave-equation-based methods, wave-equation-based methods can properly handle caustics and other problems related to ray theory. In addition, wave-equation-based methods can compute amplitudes simultaneously, thus a wave-equation-based algorithm called the SWEET (Yang et al., 2003) is used, wherein the SWEET solves the wave equation in the Laplace domain. In this algorithm, a seismic trace is considered as a series of weighted spikes. By solving the wave equation in the Laplace domain, all the spikes except the first-arrival event are attenuated and become negligible. Thus, a seismic trace originating from a series of weighted spikes can be approximated by a single spike.
As a result, the first-arrival traveltime can be extracted from the solution of wave equation in the Laplace domain. To solve the Laplace domain wave equation, a standard finite-element method is used, wherein the method exploits the combination of consistent and lumped mass matrices. The resultant matrix equation is solved by preconditioned conjugate gradient method. A numerical modeling by the Laplace domain wave equation shows less numerical dispersion error than time-domain or frequency-domain modeling (Shin et al., 2002; Shin and Cha, 2008). This characteristic enables the numerical modeling to adopt large grid spacing. However, this choice hinders exact simulation of a shallow depth source for the modeling with a coarse-grid in acoustic media. Cha and Shin (2010) used adaptive meshes to solve this problem, but implementation thereof is very complicated and the algorithm does not converge well when an iterative solver is used.
A traveltime calculation method using a conventional SWEET algorithm will be described below. The SWEET algorithm, which was suggested by Shin et al. (2002), calculates the first-arrival travel time by solving the Laplace domain wave equation. Yang et al. (2003) applied the SWEET algorithm to a 3D problem using a direct large sparse matrix solver. To calculate a first-arrival traveltime, the SWEET algorithm makes use of an assumption that a seismic signal is equal to a series of weighted spikes (see FIG. 3). Multiplying this signal by exponentially decreasing the function of time, amplitudes of all spikes except first-arrival signal are attenuated and become negligible (see FIG. 4). Therefore, the seismic signal in the time domain, which is a series of weighted spikes, can be approximated to a single spiky pulse (Shin et al., 2002). This single spiky pulse equals the first-arrival signal in the time domain. Although the principle of the SWEET algorithm can be easily explained in the time domain, the actual traveltime calculation is performed in the Laplace domain. The first-arrival traveltime is calculated using the Laplace domain wavefield and its derivative. The first-arrival traveltime t is expressed by the following equation 1:
                                          t            ⁡                          (                              x                ,                y                ,                z                ,                                  s                  opt                                            )                                =                      -                                          1                                  u                  ⁡                                      (                                          x                      ,                      y                      ,                      z                      ,                                              s                        opt                                                              )                                                              ⁡                              [                                                      ∂                                          u                      ⁡                                              (                                                  x                          ,                          y                          ,                          z                          ,                                                      s                            opt                                                                          )                                                                                                  ∂                    s                                                  ]                                                    ,                            (                  Equation          ⁢                                          ⁢          1                )            
where u is a Laplace domain wavefield, S is a Laplace domain variable, and Sopt is an optimal Laplace damping constant.
The optimal Laplace damping constant Sopt can be determined by the following equation 2 (Shin et al., 2002):
                              s          opt                =                              2            ⁢            π            ⁢                                                  ⁢                          v              ave                                            G            ⁢                                                  ⁢            Δ                                              (                  Equation          ⁢                                          ⁢          2                )            
where, νave is an average velocity of a given model, Δ is a grid spacing, and G is a number of grid points per pseudo-wavelength.
The Laplace domain wavefield is obtained by solving the 3D acoustic wave equation in the Laplace domain by the following equation 3:
                                                                                          s                  2                                                  v                  2                                            ⁢              u                        -                                                            ∂                  2                                ⁢                u                                            ∂                                  x                  2                                                      -                                                            ∂                  2                                ⁢                u                                            ∂                                  y                  2                                                      -                                                            ∂                  2                                ⁢                u                                            ∂                                  z                  2                                                              =          f                ,                            (                  Equation          ⁢                                          ⁢          3                )            
where, v is a propagation velocity in the medium, and f is a source function in the Laplace domain.
Using the finite-element method, equation 3 can be expressed as the linear algebraic system (Marfurt, 1984) by the following equations 4 and 5:Su=f  (Equation 4)withS=K+s2M,  (Equation 5)
where, S is the impedance matrix, u is the wavefield vector in the Laplace domain, f is the source vector in the Laplace domain, K is the stiffness matrix, and M is the mass matrix.
In the mean time, perfectly matched layer (PML) boundary condition is applied to eliminate unwanted edge reflections (Cohen, 2002).
To efficiently solve the equation 4, the preconditioned conjugate gradient method (Pyun et al., 2011b) is used. The partial derivative of wavefield in equation 1 is calculated by a back-propagation algorithm (Shin et al., 2002) by the following equation 6:
                                          ∂            u                                ∂            s                          =                                            S                              -                1                                      ⁡                          (                                                -                                                            ∂                      S                                                              ∂                      s                                                                      ⁢                u                            )                                .                                    (                  Equation          ⁢                                          ⁢          6                )            
The traveltime calculation method configured as above using the conventional SWEET algorithm can accurately calculate a wavefield when a source is located at a grid point with an assumption that grid spacing is small enough to avoid numerical dispersion. In particular, the Laplace domain wave equation allows accurate modeling for relatively large grid spacing when the source is located at a grid point. However, when a real source is located at shallow depth close to free surface, a problem arises that the wavefield cannot be accurately calculated by using coarse-grid spacing.
The foregoing is intended merely to aid in the understanding of the background of the present invention, and is not intended to mean that the present invention falls within the purview of the related art that is already known to those skilled in the art.